Description
Libraries
library(here)
library(data.table)
library(magrittr)
library(ggplot2)
library(scales)
library(colorspace)
library(janitor)
theme_set(theme_bw())
Making maps with R
Data from the North Carolina data
library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
# transform
d <- st_as_sf(d)
# infant deaths 1974
d$vble <- d$SID74
# infant deaths 1979
d$vble2 <- d$SID79
head(d)
## Simple feature collection with 6 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## Geodetic CRS: NAD27
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## NWBIR74 BIR79 SID79 NWBIR79 geometry vble vble2
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3... 1 0
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3... 0 3
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3... 5 6
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3... 1 2
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3... 9 3
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3... 7 5
ggplot2
library(ggplot2)
library(viridis)
# ggplot(d) +
# geom_sf(aes(fill = vble)) +
# scale_fill_viridis() +
# theme_bw()
ggplot(d) +
geom_sf(aes(fill = vble)) +
scale_fill_continuous_sequential() +
theme_bw() +
theme(
legend.position = "top"
)

plotly
library(plotly)
g <- ggplot(d) +
geom_sf(aes(fill = vble), color = "gray5", size = 0.1) +
scale_fill_continuous_sequential() +
theme_bw()
ggplotly(g)
leaflet
The sf object that we pass to leaflet() needs to have a geographic coordinate reference system (CRS) (EPSG code 4326) indicating latitude and longitude. Here, we use the st_transform() function of sf to transform the data d which has CRS given by EPSG code 4267 to CRS with EPSG code 4326.
st_crs(d)
## Coordinate Reference System:
## User input: NAD27
## wkt:
## GEOGCRS["NAD27",
## DATUM["North American Datum 1927",
## ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4267]]
library(leaflet)
pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)
l <- leaflet(d) %>%
addTiles() %>%
addPolygons(
color = "white",
fillColor = ~ pal(vble),
fillOpacity = 0.8
) %>%
addLegend(pal = pal,
values = ~ vble,
opacity = 0.8)
l
l %>% addMiniMap()
mapview
library(mapview)
mapview(d, zcol = "vble")
library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# version 2
mapview(d,
zcol = "vble",
map.types = "CartoDB.DarkMatter",
col.regions = pal)
# version 3
map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)
Side by side plots with mapview
library(leaflet.extras2) # similar to patchwork
m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m1 | m2
Synchronized maps with leafsync
m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m <- leafsync::sync(m1, m2)
m
htmltools::save_html(html = m,
file = here::here("results/figures", paste(Sys.Date(), "map-side-by-side.html", sep = "_")))
tmap
library(tmap)
# define the mode
tmap_mode("plot") # static
# tmap_mode("view") # interactive
tm_shape(d) +
tm_polygons("vble") #+

# tm_shape(st_centroid(d)) +
# tm_dots("vble")
Mobility flows with flowmapblue
- To map mobility data
- data from cellphone companies
# devtools::install_github("FlowmapBlue/flowmapblue.R")
library(flowmapblue)
# data frame 1
locations <- data.frame(
id = c(1, 2, 3),
name = c("New York", "London", "Rio de Janeiro"),
lat = c(40.713543, 51.507425, -22.906241),
lon = c(-74.011219, -0.127738, -43.180244)
)
locations
## id name lat lon
## 1 1 New York 40.71354 -74.011219
## 2 2 London 51.50742 -0.127738
## 3 3 Rio de Janeiro -22.90624 -43.180244
# atadate frame 2
flows <- data.frame(
origin = c(1, 2, 3, 2, 1, 3),
dest = c(2, 1, 1, 3, 3 , 2),
count = c(42, 51, 50, 40, 22, 42)
)
flows
## origin dest count
## 1 1 2 42
## 2 2 1 51
## 3 3 1 50
## 4 2 3 40
## 5 1 3 22
## 6 3 2 42
# final plot of transitions fromand to the 3 destinations
# flowmapblue(locations = locations, flows = flows, mapboxAccessToken = NULL, clustering = TRUE, darkMode = TRUE, animation = TRUE)
LS0tDQpvdXRwdXQ6IA0KICBib29rZG93bjo6aHRtbF9kb2N1bWVudDI6DQogICAgdGhlbWU6IGZsYXRseQ0KICAgIGNzczogInN0eWxlLmNzcyINCiAgICB0b2M6IHRydWUNCiAgICB0b2NfZGVwdGg6IDMNCiAgICB0b2NfZmxvYXQ6DQogICAgICBjb2xsYXBzZWQ6IHRydWUNCiAgICAgIHNtb290aF9zY3JvbGw6IG5vDQogICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KICAgIGNvZGVfZm9sZGluZzogc2hvdw0KdGl0bGU6ICJNYWtpbmcgbWFwcyB3aXRoIFIiDQpkYXRlOiAiU3RhcnQ6IDIyLzAzLzIwMjMiDQphdXRob3I6DQogIG5hbWU6IEpvYW8gTWFsYXRvDQogIGFmZmlsaWF0aW9uOiBJbnN0aXR1dG8gZGUgTWVkaWNpbmEgTW9sZWN1bGFyICYgSW1tdW5lLVN0YXRzDQotLS0NCg0KYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsDQogICAgICAgICAgICAgICAgICAgICAgd2FybmluZyA9IEZBTFNFLA0KICAgICAgICAgICAgICAgICAgICAgIG1lc3NhZ2UgPSBGQUxTRSwgDQogICAgICAgICAgICAgICAgICAgICAgZmlnLmFsaWduID0gImNlbnRlciIsDQogICAgICAgICAgICAgICAgICAgICAgZmlnLmFzcCA9IDAuNjE4LA0KICAgICAgICAgICAgICAgICAgICAgIGZpZy53aWR0aCA9IDEwLA0KICAgICAgICAgICAgICAgICAgICAgIGRwaSA9IDEyMCwgDQogICAgICAgICAgICAgICAgICAgICAgb3V0LndpZHRoID0gIjc1JSIpDQpgYGANCg0KIyBEZXNjcmlwdGlvbiB7LX0NCg0KDQoNCiMgTGlicmFyaWVzIHstfQ0KDQpgYGB7cn0NCmxpYnJhcnkoaGVyZSkNCmxpYnJhcnkoZGF0YS50YWJsZSkNCmxpYnJhcnkobWFncml0dHIpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KHNjYWxlcykNCmxpYnJhcnkoY29sb3JzcGFjZSkNCmxpYnJhcnkoamFuaXRvcikNCnRoZW1lX3NldCh0aGVtZV9idygpKQ0KDQpgYGANCg0KIyBNYWtpbmcgbWFwcyB3aXRoIFINCg0KLSBbbGlua10oaHR0cHM6Ly93d3cucGF1bGFtb3JhZ2EuY29tL2Jvb2stZ2RzLzI1LXNwYXRpYWwtbWFraW5nbWFwcy5odG1sKQ0KDQpEYXRhIGZyb20gdGhlIE5vcnRoIENhcm9saW5hIGRhdGENCg0KYGBge3J9DQpsaWJyYXJ5KHNmKQ0KbmFtZXNocCA8LSBzeXN0ZW0uZmlsZSgic2hhcGUvbmMuc2hwIiwgcGFja2FnZSA9ICJzZiIpDQpkIDwtIHN0X3JlYWQobmFtZXNocCwgcXVpZXQgPSBUUlVFKQ0KIyB0cmFuc2Zvcm0NCmQgPC0gc3RfYXNfc2YoZCkNCg0KIyBpbmZhbnQgZGVhdGhzIDE5NzQNCmQkdmJsZSA8LSBkJFNJRDc0DQojIGluZmFudCBkZWF0aHMgMTk3OQ0KZCR2YmxlMiA8LSBkJFNJRDc5DQpoZWFkKGQpDQpgYGANCg0KIyBnZ3Bsb3QyDQoNCmBgYHtyfQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeSh2aXJpZGlzKQ0KIyBnZ3Bsb3QoZCkgKw0KIyAgIGdlb21fc2YoYWVzKGZpbGwgPSB2YmxlKSkgKw0KIyAgIHNjYWxlX2ZpbGxfdmlyaWRpcygpICsgDQojICAgdGhlbWVfYncoKQ0KZ2dwbG90KGQpICsNCiAgZ2VvbV9zZihhZXMoZmlsbCA9IHZibGUpKSArDQogIHNjYWxlX2ZpbGxfY29udGludW91c19zZXF1ZW50aWFsKCkgKyANCiAgdGhlbWVfYncoKSArDQogIHRoZW1lKA0KICAgIGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiDQogICkNCmBgYA0KDQojIyBwbG90bHkNCg0KYGBge3J9DQpsaWJyYXJ5KHBsb3RseSkNCmcgPC0gZ2dwbG90KGQpICsNCiAgZ2VvbV9zZihhZXMoZmlsbCA9IHZibGUpLCBjb2xvciA9ICJncmF5NSIsIHNpemUgPSAwLjEpICsNCiAgc2NhbGVfZmlsbF9jb250aW51b3VzX3NlcXVlbnRpYWwoKSArIA0KICB0aGVtZV9idygpDQpnZ3Bsb3RseShnKQ0KYGBgDQoNCg0KIyBsZWFmbGV0DQoNClRoZSBzZiBvYmplY3QgdGhhdCB3ZSBwYXNzIHRvIF9sZWFmbGV0KClfIG5lZWRzIHRvIGhhdmUgYSBnZW9ncmFwaGljIGNvb3JkaW5hdGUgcmVmZXJlbmNlIHN5c3RlbSAoQ1JTKSAoRVBTRyBjb2RlIDQzMjYpIGluZGljYXRpbmcgbGF0aXR1ZGUgYW5kIGxvbmdpdHVkZS4gSGVyZSwgd2UgdXNlIHRoZSBfc3RfdHJhbnNmb3JtKClfIGZ1bmN0aW9uIG9mIGBzZmAgdG8gdHJhbnNmb3JtIHRoZSBkYXRhIGQgd2hpY2ggaGFzIENSUyBnaXZlbiBieSBFUFNHIGNvZGUgNDI2NyB0byBDUlMgd2l0aCBFUFNHIGNvZGUgNDMyNi4NCg0KDQoNCmBgYHtyfQ0Kc3RfY3JzKGQpDQpgYGANCg0KYGBge3J9DQpsaWJyYXJ5KGxlYWZsZXQpDQoNCnBhbCA8LSBjb2xvck51bWVyaWMocGFsZXR0ZSA9ICJZbE9yUmQiLCBkb21haW4gPSBkJHZibGUpDQoNCmwgPC0gbGVhZmxldChkKSAlPiUgDQogIGFkZFRpbGVzKCkgJT4lDQogIGFkZFBvbHlnb25zKA0KICAgIGNvbG9yID0gIndoaXRlIiwNCiAgICBmaWxsQ29sb3IgPSB+IHBhbCh2YmxlKSwNCiAgICBmaWxsT3BhY2l0eSA9IDAuOA0KICApICU+JQ0KICBhZGRMZWdlbmQocGFsID0gcGFsLA0KICAgICAgICAgICAgdmFsdWVzID0gfiB2YmxlLA0KICAgICAgICAgICAgb3BhY2l0eSA9IDAuOCkNCmwNCmBgYA0KDQoNCmBgYHtyfQ0KbCAlPiUgYWRkTWluaU1hcCgpDQpgYGANCg0KIyMgU2F2ZSBpbnRlcmFjdGl2ZSBmaWd1cmUgYXMgcG5nDQoNCmBgYHtyLCBldmFsPUZBTFNFfQ0KIyBTYXZlcyBtYXAuaHRtbA0KDQpsaWJyYXJ5KGh0bWx3aWRnZXRzKQ0Kc2F2ZVdpZGdldCh3aWRnZXQgPSBsLCBmaWxlID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC5odG1sIiwgc2VwID0gIl8iKSkpDQoNCiMgVGFrZXMgYSBzY3JlZW5zaG90IG9mIHRoZSBtYXAuaHRtbCBjcmVhdGVkIGFuZCBzYXZlcyBpdCBhcyBtYXAucG5nDQpsaWJyYXJ5KHdlYnNob3QpDQojIHdlYnNob3Q6Omluc3RhbGxfcGhhbnRvbWpzKCkNCndlYnNob3QodXJsID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC5odG1sIiwgc2VwID0gIl8iKSksIGZpbGUgPSBoZXJlOjpoZXJlKCJyZXN1bHRzL2ZpZ3VyZXMiLCBwYXN0ZShTeXMuRGF0ZSgpLCAibWFwLnBuZyIsIHNlcCA9ICJfIikpKQ0KYGBgDQoNCg0KIyBtYXB2aWV3DQoNCmBgYHtyfQ0KbGlicmFyeShtYXB2aWV3KQ0KbWFwdmlldyhkLCB6Y29sID0gInZibGUiKQ0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeShSQ29sb3JCcmV3ZXIpDQpwYWwgPC0gY29sb3JSYW1wUGFsZXR0ZShicmV3ZXIucGFsKDksICJZbE9yUmQiKSkNCiMgdmVyc2lvbiAyDQptYXB2aWV3KGQsDQogICAgICAgIHpjb2wgPSAidmJsZSIsDQogICAgICAgIG1hcC50eXBlcyA9ICJDYXJ0b0RCLkRhcmtNYXR0ZXIiLA0KICAgICAgICBjb2wucmVnaW9ucyA9IHBhbCkNCmBgYA0KDQpgYGB7cn0NCiMgdmVyc2lvbiAzDQptYXAxIDwtIG1hcHZpZXcoZCwgemNvbCA9ICJ2YmxlIikNCmxlYWZsZXQ6OmFkZE1pbmlNYXAobWFwMUBtYXApDQpgYGANCg0KDQoNCiMgU2lkZSBieSBzaWRlIHBsb3RzIHdpdGggbWFwdmlldw0KDQpgYGB7cn0NCmxpYnJhcnkobGVhZmxldC5leHRyYXMyKSAjIHNpbWlsYXIgdG8gcGF0Y2h3b3JrDQptMSA8LSBtYXB2aWV3KGQsIHpjb2wgPSAidmJsZSIpDQptMiA8LSBtYXB2aWV3KGQsIHpjb2wgPSAidmJsZTIiKQ0KbTEgfCBtMg0KYGBgDQoNCg0KIyBTeW5jaHJvbml6ZWQgbWFwcyB3aXRoIGxlYWZzeW5jDQoNCmBgYHtyfQ0KbTEgPC0gbWFwdmlldyhkLCB6Y29sID0gInZibGUiKQ0KbTIgPC0gbWFwdmlldyhkLCB6Y29sID0gInZibGUyIikNCm0gPC0gbGVhZnN5bmM6OnN5bmMobTEsIG0yKQ0KbQ0KYGBgDQoNCg0KYGBge3IsIGV2YWw9RkFMU0V9DQpodG1sdG9vbHM6OnNhdmVfaHRtbChodG1sID0gbSwgDQogICAgICAgICAgICAgICAgICAgICBmaWxlID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC1zaWRlLWJ5LXNpZGUuaHRtbCIsIHNlcCA9ICJfIikpKQ0KYGBgDQoNCg0KIyBgdG1hcGANCg0KYGBge3J9DQpsaWJyYXJ5KHRtYXApDQoNCiMgZGVmaW5lIHRoZSBtb2RlDQp0bWFwX21vZGUoInBsb3QiKSAjIHN0YXRpYw0KIyB0bWFwX21vZGUoInZpZXciKSAjIGludGVyYWN0aXZlDQoNCnRtX3NoYXBlKGQpICsgDQogIHRtX3BvbHlnb25zKCJ2YmxlIikgIysNCiAgIyB0bV9zaGFwZShzdF9jZW50cm9pZChkKSkgKyANCiAgIyB0bV9kb3RzKCJ2YmxlIikNCmBgYA0KDQojIE1vYmlsaXR5IGZsb3dzIHdpdGggYGZsb3dtYXBibHVlYA0KDQorIFRvIG1hcCBtb2JpbGl0eSBkYXRhDQorIGRhdGEgZnJvbSBjZWxscGhvbmUgY29tcGFuaWVzDQoNCmBgYHtyfQ0KIyBkZXZ0b29sczo6aW5zdGFsbF9naXRodWIoIkZsb3dtYXBCbHVlL2Zsb3dtYXBibHVlLlIiKQ0KbGlicmFyeShmbG93bWFwYmx1ZSkNCmBgYA0KDQpgYGB7cn0NCiMgZGF0YSBmcmFtZSAxDQpsb2NhdGlvbnMgPC0gZGF0YS5mcmFtZSgNCiAgaWQgPSBjKDEsIDIsIDMpLA0KICBuYW1lID0gYygiTmV3IFlvcmsiLCAiTG9uZG9uIiwgIlJpbyBkZSBKYW5laXJvIiksDQogIGxhdCA9IGMoNDAuNzEzNTQzLCA1MS41MDc0MjUsIC0yMi45MDYyNDEpLA0KICBsb24gPSBjKC03NC4wMTEyMTksIC0wLjEyNzczOCwgLTQzLjE4MDI0NCkNCikNCmxvY2F0aW9ucw0KDQojIGF0YWRhdGUgZnJhbWUgMg0KZmxvd3MgPC0gZGF0YS5mcmFtZSgNCiAgb3JpZ2luID0gYygxLCAyLCAzLCAyLCAxLCAzKSwNCiAgZGVzdCA9IGMoMiwgMSwgMSwgMywgMyAsIDIpLA0KICBjb3VudCA9IGMoNDIsIDUxLCA1MCwgNDAsIDIyLCA0MikNCikNCmZsb3dzDQpgYGANCg0KYGBge3J9DQojIGZpbmFsIHBsb3Qgb2YgdHJhbnNpdGlvbnMgZnJvbWFuZCB0byB0aGUgMyBkZXN0aW5hdGlvbnMNCiMgZmxvd21hcGJsdWUobG9jYXRpb25zID0gbG9jYXRpb25zLCBmbG93cyA9IGZsb3dzLCBtYXBib3hBY2Nlc3NUb2tlbiA9IE5VTEwsIGNsdXN0ZXJpbmcgPSBUUlVFLCBkYXJrTW9kZSA9IFRSVUUsIGFuaW1hdGlvbiA9IFRSVUUpDQpgYGANCg0KDQoNCg==